function out = iterative_soft_thresholding( varargin )
% out = iterative_soft_thresholding[Y,angles,truncate,imgSize,maxit,intialGuessStr,constThreshValue]
%
% Juergen Frikel 25/11/2011
%
% Sparse recovery: Accelerated iterated soft thresholding (AIST) with
%                  steplength selection and adapted universal threshold
%
% INPUT:
%
%   Y                   [noisy] projection Data 
%   
%   angles              measurement angles
%   
%   truncate            (optional) boolean flag. If its value is 'true', the dimensionality
%                       reduction in the curvelet domain is performed.
%
%   imgSize             (optional) size of the image to be reconstructed
%
%   maxit               (optional) max. number of iteration
%
%   initialGuessStr     (optional) initial guess
%
%   constThreshValue    regularization (thresholding) parameter
%
%
% OUTPUT variable 'out' is a struct array containing the following entries
%
%   rec                 reconstructed image
%
%   visibleCurveletIdx  indices of visible curvelets at the angular range
%                       given by 'angles'
%
%   iterates            cell array containing all iterates
%
%   error               vector of data errors of the corresponding iterates
%      
%
% Literature:
%   I. Loris, M. Bertero, C. De Mol, R. Zanella, L. Zanni. "Accelerating 
%   gradient projection methods for 1-constrained signal recovery by 
%   steplength selection rules", Appl. Comput. Harmon. Anal. 27 (2009),
%   247-254.


% fprintf(1,'Main Routine: Accelerated iterated soft thresholding (AIST)\n\n');
% fprintf(1,'\nProjection of CURVELET coefficients\n\n');

% --------------- LOOKUP TABLE for l^p-thresholding -----------------------
global x_thresh y_thresh p_thresh

x_thresh = linspace(-10,10,1000);
y_thresh = [];
p_thresh = 1.1;
% -------------------------------------------------------------------------

% intialize output variable
out = struct( 'rec',[],'visibleCurveletIdx',[],'iterates',[],'error',[],'rectime',[] );

% Main computation
% [coeffs,ccindex,err,iterates,N] = iterate_with_steplengthSelectionForGP( varargin{:} );
tic;

[coeffs,out.visibleCurveletIdx,out.error,out.iterates,N] = iterative_softthresholding( varargin{:} );

out.rec     = ifdct( coeffs, out.visibleCurveletIdx, N, N);
out.rectime = toc;

% sec = round(toc);
% s   = mod(sec,60);
% m   = (sec-s)/60;
% fprintf(1,'Elapsed Time: %2d min %2d sec\n\n',m,s);

%==========================================================================
function y = prox_lp( x, weight, p )

% % given x, y is computed by solving F(x,y)=0 via Newton method
% F   = @(x,y)(y^4-2*x*y^3+x^2*y^2+epsilon^2*y^2-weight^2*y^2-2*epsilon^2*x*y+epsilon^2*x^2);
% dyF = @(x,y)(4*y^3-6*x*y^2+2*x^2*y+2*epsilon^2*y-2*weight^2*y-2*epsilon^2*x);
% 
% % choose intial guess as y0 = x;
% maxIt = 100;
% nofIt = 0;
% y = x;
% while( F(x,y)>1e-8 || nofIt<=maxIt )
%    y = y - F(x,y)/dyF(x,y);
%    
%    nofIt = nofIt + 1;
% end

F = @(y)(x-y-weight*p*sign(y)*abs(y)^(p-1));
y = fsolve(F,x);

%==========================================================================
function z = addCurveletCoeffs( x, factor, y, ccindex )
% compute: z = x + factor*y

lx = length(x);
ly = length(y);

if( lx~=ly )
    error('Coeff lengths must be equal!');
end

% allocate memory for the resulting coeff array
z = x;

% add approximation coeffs
if( size(x{1}{1})~=size(y{1}{1}) )
    error('Coeffs must have equal coeff sizes');
else
    z{1}{1} = x{1}{1} + factor.*y{1}{1};
end

for k=2:lx
   llx = length( x{k} );
   lly = length( y{k} );
   
   if( llx~=lly )
       error('Coeffs must have equal sizes');
   end
   
   for n=ccindex{k-1}
      if( size(x{k}{n})~=size(y{k}{n}) )
          error('Coeffs must have equal coeff sizes');
      end
      
      z{k}{n} = x{k}{n} + factor.*y{k}{n};
   end
end;
% z{end}{1} = x{end}{1} + factor.*y{end}{1};
%==========================================================================

function sigma = estimate_std( coeffs, ccindex )

tmp = [];
for n=ccindex{end}
    tmp = [tmp; abs(coeffs{end}{n}(:))];
end

sigma = median(tmp)/0.6745;
fprintf(1,'\nfinest scale coefficients: min = %e\tmax = %e\tmedian = %e\tsigma = %e\n',min(tmp),max(tmp),median(tmp),sigma);

%==========================================================================

function [step,thresh] = steplegth( x, ccindex, dimRed )

sigma  = (1-dimRed)*estimate_std( x, ccindex );
maxscale = length(x)+1;

elPerScale = zeros(1,length(x)-1);
thresh = zeros(1,length(x)-1);
for j=2:length(x)
    for l=1:length(x{j})
       elPerScale(j-1) = elPerScale(j-1) + numel(x{j}{l});
    end
    
    thresh(j-1) = 2^((j-maxscale)*3/4)*sigma*sqrt(2*log(elPerScale(j-1)));
end

step = 2*min( thresh );
fprintf(1,'sigma = %e\tsteplenght = %e\n\n',step,sigma);

%==========================================================================

function [x,l1norm] = softThresh_curvelet2( x, ccindex, threshvalue )

l1norm = 0;
for k=2:length(x)
    for n=ccindex{k-1}
        x{k}{n} = x{k}{n}.* (abs(x{k}{n}) > threshvalue(k-1));
        x{k}{n} = x{k}{n} - sign(x{k}{n}).*threshvalue(k-1);
        
        l1norm = l1norm + sum( sum( abs(x{k}{n}) ) );
    end
end

%==========================================================================

function threshSeq = getThresholdingSequence( x, ccindex, dimRed )

sigma  = (1-dimRed)*estimate_std( x, ccindex );

if( sigma==0 )
    disp('NO TRESHOLDING!!!');
end

l1norm = 0;
maxscale = length(x)+1;
threshSeq = cell(size(x));
for k=2:length(x)
    %thres = 2^((k-maxscale)*3/4)*sigma;
    thres = 2^((k-maxscale))*sigma;
    for n=ccindex{k-1}
        lambda = thres*sqrt(2*log(numel(x{k}{n})));
        %lambda = thres^2*sqrt(2*log(numel(x{k}{n})));
        
        threshSeq{k}{n} = lambda;

        x{k}{n} = x{k}{n}.* (abs(x{k}{n}) > lambda);
        x{k}{n} = x{k}{n} - sign(x{k}{n}).*lambda;
        
        l1norm = l1norm + sum( sum( abs(x{k}{n}) ) );
    end
end

%==========================================================================

function [x,l1norm] = softThresh_curvelet( x, ccindex, threshSeq )

l1norm = 0;
for k=2:length(x)
    for n=ccindex{k-1}
        x{k}{n} = x{k}{n}.* (abs(x{k}{n}) > threshSeq{k}{n});
        x{k}{n} = x{k}{n} - sign(x{k}{n}).* threshSeq{k}{n};
        
        l1norm = l1norm + sum( sum( abs(x{k}{n}) ) );
    end
end

%==========================================================================

function [x,l1norm] = softThresh_curvelet_const( x, ccindex, lambda )

l1norm = 0;
for k=2:length(x)
    for n=ccindex{k-1}
        x{k}{n} = x{k}{n}.* (abs(x{k}{n}) > lambda);
        x{k}{n} = x{k}{n} - sign(x{k}{n}).* lambda;
        
        l1norm = l1norm + sum( sum( abs(x{k}{n}) ) );
    end
end


%==========================================================================

function [x,l1norm] = generic_softThresh_curvelet( x, ccindex, dimRed )

sigma  = (1-dimRed)*estimate_std( x, ccindex );

if( sigma==0 )
    disp('NO TRESHOLDING!!!');
end

l1norm = 0;
maxscale = length(x)+1;
for k=2:length(x)
    %thres = 2^((k-maxscale)*3/4)*sigma;
    thres = 2^((k-maxscale))*sigma;
    for n=ccindex{k-1}
        lambda = thres*sqrt(2*log(numel(x{k}{n})));
        %lambda = thres^2*sqrt(2*log(numel(x{k}{n})));

        x{k}{n} = x{k}{n}.* (abs(x{k}{n}) > lambda);
        x{k}{n} = x{k}{n} - sign(x{k}{n}).*lambda;
        
        l1norm = l1norm + sum( sum( abs(x{k}{n}) ) );
    end
end

%==========================================================================
function vector = extractCoeffs( coeffs, ccindex )
% extract visible curvelet coeffs to vector

% approximation coeffs
vector = coeffs{1}{1}(:);

% visible detail coeffs
for k=2:length(coeffs)
    for n=ccindex{k-1}
        vector = [vector;coeffs{k}{n}(:)];
    end
end

%vector = [vector;coeffs{end}{1}(:)];

%==========================================================================
function index = findVisibleCurveltIndices( maxscale, theta_min, theta_max )
% Input: 
%   maxscale    the finest (largest) scale for detail coeffs
%   theta_min
%   theta_max   angular range for the coeffs to be kept
%
% Output:
%   index       index array of DETAIL curvelet coeffs that should be kept
%               organization of the index array: 
%               index{1}    indices for coeffs{2}
%               index{2}    indices for coeffs{3}
%               ...
%               index{end}  indices for coeffs{end-1}

if( theta_max >= 180 )
    theta_max = 179;
end
theta_min = mod( theta_min, 180 );
theta_max = mod( theta_max, 180 );

theta_min = pi*theta_min/180;
theta_max = pi*theta_max/180;

if( maxscale < 3 )
    error('maxscale should be greater or equal than 3');
end

index = [];

tan_min = tan( theta_min );
tan_max = tan( theta_max );

if( tan_min>1e+10 )
    tan_min = Inf;
end
if( tan_max>1e+10 )
    tan_max = Inf;
end

itan_min = 1./tan_min;
itan_max = 1./tan_max;

for scale=3:maxscale

    N = 2*2^floor(scale/2);

    %fprintf(1,'tan_min = %f\n', tan_min);
    %fprintf(1,'tan_max = %f\n', tan_max);

    idx = [];
    
    idx_min = 0;
    idx_max = 0;

    dSlope = 2/N;
    % find idx_min
    if( tan_min <=1 && tan_min > -1 )
        for k=1:N
            slope1 = 1 - (k-1)*dSlope;
            slope2 = 1 - k*dSlope;

            if( tan_min<=slope1 && tan_min>slope2 )
                idx_min = k + N;
                break;
            end
        end
    else
        for k=1:N
           slope1 = -1 + (k-1)*dSlope;
           slope2 = -1 + k*dSlope;

           if( itan_min>=slope1 && itan_min<slope2 )
               idx_min = k;
               break;
           end
        end
    end
    % find idx_max
    if( tan_max <=1 && tan_max > -1 )
        for k=1:N
            slope1 = 1 - (k-1)*dSlope;
            slope2 = 1 - k*dSlope;

            if( tan_max<=slope1 && tan_max>slope2 )
                idx_max = k + N;
                break;
            end
        end
    else
        for k=1:N
           slope1 = -1 + (k-1)*dSlope;
           slope2 = -1 + k*dSlope;

           if( itan_max>=slope1 && itan_max<slope2 )
               idx_max = k;
               break;
           end
        end
    end

    % since the numbering of the curvelet coeffs starts at 135 degree, there
    % are 2 distinct regions of angles which have to be treated separately,
    % nameley [0 135] and [135 180] degrees.

    % chek if theta_min and theta_max are both within the same region
    border = 0.75*pi;
    sameregion = ( theta_min<=border && theta_max<=border ) ||...
                 ( theta_min>=border && theta_max>=border );

    tmp = idx_min;
    idx_min = min( [idx_min idx_max] );
    idx_max = max( [tmp idx_max] );
    
    % find indices of invisible curvelet-coefficients
    % Note: theta_min and theta_max describe the range of visible
    % curvelet-coeffs. Therefor, we are interested in the complement of the
    % index-set described by these angles
    if( theta_min<=theta_max)
        if( sameregion )
            % visible range is: [idx_min:idx_max]
            idx = idx_min:idx_max;
        else
            % visible range is: [1:idx_min idx_max:end]
            idx = [1:idx_min idx_max:2*N];
        end
    else
        if( sameregion )
            % visible range is: [1:idx_min idx_max:end]
            idx = [1:idx_min idx_max:2*N];
        else
            % visible range is: [idx_min:idx_max]
            idx = idx_min:idx_max;
        end
    end
    idx = mod( idx-1, 2*N ) + 1;

    
    % fill by symmetry
    index{scale-2} = [ idx idx+2*N ]; % index fuer Curvelet-Koeffs
end

index{end+1} = index{end};
%==========================================================================

function n = nofCurveletCoeffs( C, ccindex )

n = 0;
for k=1:length(C)
    if( k==1 )
        n = n + numel(C{1}{1});
    else
        for l=ccindex{k-1};%1:length(C{k})
            n = n + numel(C{k}{l});
        end
    end
end

%==========================================================================

function deg = degreeOfDimReduction( M,N,ccindex )

C = fdct_wrapping( ones(M,N) );
n = 0;
N = 0;
for k=1:length(C)
    if( k==1 )
        n = n + numel(C{1}{1});
        N = n;
    else
        for l=1:length(C{k})
            tmp = numel(C{k}{l});
            N = N + tmp;
            if( ismember(l,ccindex{k-1}) )
                n = n + tmp;
            end
        end
    end
end
deg = 1 - n/N;
deg = round(100*deg)/100;
%==========================================================================

%function [B,err,iterates] = iterate_with_steplengthSelectionForGP( Y,angles,tol,maxit,varargin)
function [wNew,ccindex,err,iterates,N] = iterate_with_steplengthSelectionForGP( varargin )

global exit_flag

% parse inputs
[Y,angles,truncate,imgSize,maxit,intialGuessStr,constThreshValue] = parseinputs( varargin{:} );

% global Parameters for iteration
tol   = 0.01;   % tolerance

exit_flag    = false; % exit main routine on user 
plotIterates = true; % plot each iterate?
plotError    = false;
%truncate     = true; % adapted curvelet transform?

% store all iterates in a 3D-array
iterates = [];%zeros(M,N,maxit);

N = imgSize;
if( isempty(N) )
    initialGuess = iradon(Y,angles);
    N = max(size(initialGuess));
end


switch lower(intialGuessStr)
    case 'zero'
        initialGuess = zeros(N);
    case 'fbp'
        initialGuess = iradon(Y,angles,N);
    case 'backprojection'
        initialGuess = backproject(Y,angles,N);
    otherwise
        error('Unknown option for the initial guess!');
end

nofPr = size(Y,1);

tv_weight = 0.001*(max(initialGuess(:))-min(initialGuess(:)));

% compute indices of visible coeffs
ccindex = [];
nbscales = ceil(log2(N) - 3);
if( truncate )
    %theta_min = 0;
    %theta_max = 90;
    
    % find curvlet indices corresponding to the angular range
    % [theta_min, theta_max]
    ccindex = findVisibleCurveltIndices( nbscales , min(angles(:)), max(angles(:)) );
    
    fprintf(1,'TRUNCATION of invisible curvelet coefficients\n');
    fprintf(1,'Visible angular range: [ %d, %d ]\n',min(angles(:)), max(angles(:)));
    fprintf(1,'Nof DETAILS = %d\n',length(cell2mat(ccindex)));
else
    % determine the full index array of all curvelet coefficients
    ccindex = findVisibleCurveltIndices( nbscales, 0, 179 );
    
    fprintf(1,'Computation of all coefficients\n'); 
    fprintf(1,'Nof DETAILS = %d\n',length(cell2mat(ccindex)));
end

% Initializations
wOld = fdct( initialGuess, ccindex );
wNew = wOld;
h    = wOld;

dimRed = degreeOfDimReduction(N,N,ccindex);

fprintf(1,'Dimension in the curvelet domain: %d\n',nofCurveletCoeffs(wOld,ccindex));
fprintf(1,'Degree of DimReduction in the curvelet domain: %d\n',dimRed);

% compute residual
img = ifdct( wOld, ccindex, N, N );
proj = radon(img,angles,nofPr);
res  = norm( Y(:) - proj(:) )^2;

% Parameters of GPSS
a_min = 10e-10;
a_max = 10e+10;
beta  = 10e-4;
gamma = 0.5;    % = theta
nofEl = 10;     % = M
tau   = 0.5;
M_alpha = 2;
alpha_BB2_vec = inf*ones(1,M_alpha+1);

% vector of the last residuals
resvec = zeros( 1, nofEl );

% vector containing the data error of each iteration
err = [];

% return value
B = [];

% descent direction of current iteration
descent = wOld;
lambda  = 1.0;
%lambda_old = lambda;

% choose initial steplength
proj = radon(img,angles,nofPr);
tmp = fdct( backproject(Y-proj,angles,N), ccindex );
tmp = addCurveletCoeffs( wOld, 1.0, tmp, ccindex );
alpha = 1/norm( extractCoeffs(tmp,ccindex), inf );
alpha = min( [alpha a_max] );
alpha = max( [a_min alpha] );
    
it = 0;
if( plotIterates )
    figHandle = figure('KeyPressFcn',@keylistener);
    set(figHandle,'Position',[50 200 500 500])
    img = real( ifdct( wOld, ccindex, N, N ) );
    imagesc( img );
    colormap gray;
    title(strcat('Iteration ',num2str(it)));
    drawnow;
    
    if( plotError )
        plotHandle = figure;
        set(plotHandle,'Position',[600 200 500 500])
    end
end
while( res>=tol && it<maxit )
    it = it + 1;
    
    %figure,imshow(real( ifdct( wOld, ccindex, N, N ) )); title('1111');
    
    %TV step
%     img  = real( ifdct( wOld, ccindex, N, N ) );
%     tv_weight = tv_weight/it;
%     img  = tvexact( img, tv_weight );
%     wOld = fdct( img, ccindex );
    
    %figure,imshow(real( ifdct( wOld, ccindex, N, N ) )); title('2222');
    %return;
    
%     negidx = find( xOld<0 );
%     xOld(negidx) = 0;

    img = ifdct( wOld, ccindex, N, N );
    proj = radon(img,angles,nofPr);
    res  = norm( Y(:) - proj(:) )^2;
    
    err( it ) = sqrt(res);
    resvec(mod(it,nofEl)+1) = res;
    
    % compute descent direction
    img  = ifdct( wOld, ccindex, N, N );
    proj = radon(img,angles,nofPr);
    bp   = backproject( Y-proj, angles, N);
    grad = fdct( -bp, ccindex );% = -R'(y-R*img)
    
    %fprintf(1,'\nnorm(wOld) = %f;  \n',norm( extractCoeffs( wOld, ccindex ) ));
    %fprintf(1,'norm(alpha) = %f;  \n', alpha);
    %fprintf(1,'norm(grad) = %f;  \n',norm( extractCoeffs( grad, ccindex ) ));
    
    wNew = addCurveletCoeffs( wOld, -alpha, grad, ccindex );% h = wOld(:) - alpha*gradOld;
    
    % Soft-thresholding: h = S( xOld -alpha*grad(xOld) );
    %[wNew,l1norm] = softThresh_curvelet( wNew, ccindex, dimRed );
    l1norm = 0;
    threshvalue = 1;
    for k=2:length(wOld)
        for n=ccindex{k-1}
            wOld{k}{n} = wOld{k}{n}.* (abs(wOld{k}{n}) > threshvalue);
            wOld{k}{n} = wOld{k}{n} - sign(wOld{k}{n}).*threshvalue;

            l1norm = l1norm + sum( sum( abs(wOld{k}{n}) ) );
        end
    end
    
%     f_max  = max( resvec );
%     lambda = 1.0;%300.0??????
%     
%     img  = ifdct( wNew, ccindex, N, N );
%     proj = radon(img,angles,nofPr);
%     res  = norm( Y(:) - proj(:) )^2;
%     
%     vgrad  = extractCoeffs( grad, ccindex );
%     tmp    = extractCoeffs( descent, ccindex );
%     uBound = f_max + beta*lambda*(vgrad'*tmp(:));
    
%     subit = 0;
%     while( res>uBound && subit<50 )
%         subit = subit + 1;
%         
%         wNew = addCurveletCoeffs( wOld, -gamma*alpha, grad, ccindex );
%         
%         img  = ifdct( wNew, ccindex, N, N );
%         
%         proj = radon(img,angles,nofPr);
%         res  = norm( Y(:) - proj(:) )^2;
%         
%         vgrad = extractCoeffs( grad, ccindex );
%         tmp  = extractCoeffs( descent, ccindex );
%         uBound = f_max + beta*lambda*(vgrad'*tmp(:));
%     end

%%%%%%%%%%%%%%%%%%%%%%%%  STEPLENGTH SELECTION  %%%%%%%%%%%%%%%%%%%%%%%%%%%   
    %fprintf(1,'alpha  = %f \n',alpha);
    %fprintf(1,'lambda = %f \n\n',lambda);
        
    % Compute the 'alpha' for the next iteration
    dx  = addCurveletCoeffs( wNew, -1.0, wOld, ccindex );
    img = ifdct( dx, ccindex, N, N );
    %figure,imagesc(real(img)),colormap gray, colorbar
    dg  = radon(img,angles,nofPr);
    
    prod_dxdg = dg(:)'*dg(:);
    %fprintf(1,'dx*dg    = %e \n',prod_dxdg);
    
    if( prod_dxdg<=0 )
        alpha_new = a_max;
    else
        dg = backproject(dg,angles,N);

        dx = extractCoeffs( dx, ccindex );
        
        a_BB1 = min( [(dx(:)'*dx(:)/prod_dxdg) a_max] );
        a_BB1 = max( [a_min a_BB1] );

        a_BB2 = min( [(prod_dxdg/(dg(:)'*dg(:))) a_max] );
        a_BB2 = max( [a_min a_BB2] );
        
        %fprintf(1,'norm(dx) = %f \n',norm(dx(:))^2);
        %fprintf(1,'norm(dg) = %f \n\n',norm(dg(:))^2);
        
        %fprintf(1,'a_BB1 = %f \n',a_BB1);
        %fprintf(1,'a_BB2 = %f \n\n',a_BB2);

        % store the a_BB2
        alpha_BB2_vec( mod(it,M_alpha+1)+1 ) = a_BB2;

        if( a_BB2/a_BB1 <= tau )
            %disp('BB2')
            alpha_new = min( alpha_BB2_vec );
            tau = 0.9*tau;
        else
            %disp('BB1')
            alpha_new = a_BB1;
            tau = tau*1.1;
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   
    % new function value
    %l1norm = 0;%normCurvletCoeffs(wOld,1);
    fval = res + l1norm;
    
    if(mod(it,1)==0)
        if( plotIterates )
            %iterate = real( ifdct( wNew, ccindex, N, N ) );
            figure(figHandle);
            img = real( ifdct( wNew, ccindex, N, N ) );
            imagesc( img );
            colormap gray;
            title(strcat('Iteration ',num2str(it)));
            drawnow;
            %pause(1)
            
            if( plotError )
                figure(plotHandle)
                semilogy(it,res,'x');
                xlabel('iterations');
                ylabel('absolute error');
                hold on;
            end
        end
        fprintf(1,'it=%5d;\terr=%e;\tfval=%e;\tl1-norm=%e;\tstep=%e;\ttv_weight = %e;\n',it,err(it),fval,l1norm,alpha,tv_weight);
    end
    
%     if( mod(it,20)==0 )
%         save( strcat('iterates',num2str(it),'.mat'),'iterates'); 
%     end
    
    wOld = wNew;
    alpha = alpha_new;
    %lambda_old = lambda;
    
    % save current iterate
    iterates(:,:,it) = real( ifdct( wOld, ccindex, N, N ) )';
    
    if( exit_flag )
        return;
    end
end

%B = real( ifdct( wNew, ccindex, N, N ) );

%fprintf(1,'it=%5d;\terr=%e;\tfval=%e;\tl1-norm=%e;\tstep=%e;\tdiff=%e;\n',it,err(it),fval,l1norm,alpha,diff);
fprintf(1,'Eof Main Computation\n\n');
%Eof

% =========================================================================
function [wNew,ccindex,err,iterates,N] = iterative_softthresholding( varargin )

global exit_flag

% parse inputs
[Y,angles,truncate,imgSize,maxit,intialGuessStr,constThreshValue] = parseinputs( varargin{:} );

% global Parameters for iteration
tol   = 0.01;   % tolerance

exit_flag    = false; % exit main routine on user 
plotIterates = true; % plot each iterate?
plotError    = false;
printLog     = true;
genericThreshold = false;
constThreshold = false;

% store all iterates in a 3D-array
iterates = [];%zeros(M,N,maxit);

N = imgSize;
if( isempty(N) )
    initialGuess = iradon(Y,angles);
    N = max(size(initialGuess));
end


switch lower(intialGuessStr)
    case 'zero'
        initialGuess = zeros(N);
    case 'fbp'
        initialGuess = iradon(Y,angles,N);
    case 'backprojection'
        initialGuess = backproject(Y,angles,N);
    otherwise
        error('Unknown option for the initial guess!');
end

nofPr = size(Y,1);

tv_weight = 0.001*(max(initialGuess(:))-min(initialGuess(:)));

% compute indices of visible coeffs
ccindex = [];
nbscales = ceil(log2(N) - 3);
if( truncate )
    %theta_min = 0;
    %theta_max = 90;
    
    % find curvlet indices corresponding to the angular range
    % [theta_min, theta_max]
    ccindex = findVisibleCurveltIndices( nbscales , min(angles(:)), max(angles(:)) );
    
    if( printLog )
        fprintf(1,'TRUNCATION of invisible curvelet coefficients\n');
        fprintf(1,'Visible angular range: [ %d, %d ]\n',min(angles(:)), max(angles(:)));
        fprintf(1,'Nof DETAILS = %d\n',length(cell2mat(ccindex)));
        fprintf(1,'MaxIterations = %d\n',maxit);
        fprintf(1,'Initial guess = %s\n',intialGuessStr);
    end
else
    % determine the full index array of all curvelet coefficients
    ccindex = findVisibleCurveltIndices( nbscales, 0, 179 );
    
    if( printLog )
        fprintf(1,'Computation of all coefficients\n'); 
        fprintf(1,'Nof DETAILS = %d\n',length(cell2mat(ccindex)));
    end
end

% Initializations
wOld = fdct( initialGuess, ccindex );
wNew = wOld;
h    = wOld;

dimRed = degreeOfDimReduction(N,N,ccindex);

threshSeq = [];
if( ~genericThreshold && ~constThreshold)
    threshSeq = getThresholdingSequence( wNew, ccindex, dimRed );
end
% if( genericThreshold )
%     [wNew,l1norm] = generic_softThresh_curvelet( wNew, ccindex, 0 );
% elseif( constThreshold )
%     [wNew,l1norm] = softThresh_curvelet_const( wNew, ccindex, constThreshValue );
% else
%     [wNew,l1norm] = softThresh_curvelet( wNew, ccindex, threshSeq );
% end


if( printLog )
    fprintf(1,'Dimension in the curvelet domain: %d\n',nofCurveletCoeffs(wOld,ccindex));
    fprintf(1,'Degree of DimReduction in the curvelet domain: %d\n',dimRed);
end

% compute residual
img = ifdct( wOld, ccindex, N, N );
proj = radon(img,angles,nofPr);
res  = norm( Y(:) - proj(:) )^2;

% Parameters of GPSS
a_min = 10e-10;
a_max = 10e+10;
beta  = 10e-4;
gamma = 0.5;    % = theta
nofEl = 10;     % = M
tau   = 0.5;
M_alpha = 2;
alpha_BB2_vec = inf*ones(1,M_alpha+1);

% vector of the last residuals
resvec = zeros( 1, nofEl );

% vector containing the data error of each iteration
err = [];

% return value
B = [];

% descent direction of current iteration
descent = wOld;
lambda  = 1.0;
%lambda_old = lambda;

% choose initial steplength
proj = radon(img,angles,nofPr);
tmp = fdct( backproject(Y-proj,angles,N), ccindex );
tmp = addCurveletCoeffs( wOld, 1.0, tmp, ccindex );
alpha = 1/norm( extractCoeffs(tmp,ccindex), inf );
alpha = min( [alpha a_max] );
alpha = max( [a_min alpha] );
    
it = 0;
if( plotIterates )
    figHandle = figure('KeyPressFcn',@keylistener);
    set(figHandle,'Position',[50 200 500 500])
    img = real( ifdct( wOld, ccindex, N, N ) );
    imagesc( img );
    colormap gray;
    title(strcat('Iteration ',num2str(it)));
    drawnow;
    
    if( plotError )
        plotHandle = figure;
        set(plotHandle,'Position',[600 200 500 500])
    end
end
while( res>=tol && it<maxit )
    it = it + 1;
    
    %figure,imshow(real( ifdct( wOld, ccindex, N, N ) )); title('1111');
    
    %TV step
%     img  = real( ifdct( wOld, ccindex, N, N ) );
%     tv_weight = tv_weight/it;
%     img  = tvexact( img, tv_weight );
%     wOld = fdct( img, ccindex );
    
    %figure,imshow(real( ifdct( wOld, ccindex, N, N ) )); title('2222');
    %return;
    
%     negidx = find( xOld<0 );
%     xOld(negidx) = 0;

    img = ifdct( wOld, ccindex, N, N );
    proj = radon(img,angles,nofPr);
    res  = norm( Y(:) - proj(:) )^2;
    
    err( it ) = sqrt(res);
    resvec(mod(it,nofEl)+1) = res;
    
    % compute descent direction
%     img  = ifdct( wOld, ccindex, N, N );
%     proj = radon(img,angles,nofPr);
    bp   = backproject( Y-proj, angles, N);
    grad = fdct( -bp, ccindex );% = -R'(y-R*img)
    
    %fprintf(1,'\nnorm(wOld) = %f;  \n',norm( extractCoeffs( wOld, ccindex ) ));
    %fprintf(1,'norm(alpha) = %f;  \n', alpha);
    %fprintf(1,'norm(grad) = %f;  \n',norm( extractCoeffs( grad, ccindex ) ));
    
    wNew = addCurveletCoeffs( wOld, -alpha, grad, ccindex );% h = wOld(:) - alpha*gradOld;
    
    % Soft-thresholding: h = S( xOld -alpha*grad(xOld) );
    if( genericThreshold )
        [wNew,l1norm] = generic_softThresh_curvelet( wNew, ccindex, 0 );
    elseif( constThreshold )
        [wNew,l1norm] = softThresh_curvelet_const( wNew, ccindex, constThreshValue );
    else
        [wNew,l1norm] = softThresh_curvelet( wNew, ccindex, threshSeq );
    end
       
    % new function value
    %l1norm = 0;%normCurvletCoeffs(wOld,1);
    fval = res + l1norm;
    
    if(mod(it,1)==0)
        if( plotIterates )
            %iterate = real( ifdct( wNew, ccindex, N, N ) );
            figure(figHandle);
            img = real( ifdct( wNew, ccindex, N, N ) );
            imagesc( img );
            colormap gray;
            title(strcat('Iteration ',num2str(it)));
            drawnow;
            %pause(1)
            
            if( plotError )
                figure(plotHandle)
                semilogy(it,res,'x');
                xlabel('iterations');
                ylabel('absolute error');
                hold on;
            end
        end
        if( printLog )
            fprintf(1,'it=%5d;\terr=%e;\tfval=%e;\tl1-norm=%e;\tstep=%e;\ttv_weight = %e;\n',it,err(it),fval,l1norm,alpha,tv_weight);
        end
    end
    
%     if( mod(it,20)==0 )
%         save( strcat('iterates',num2str(it),'.mat'),'iterates'); 
%     end
    
    wOld = wNew;
    %alpha = alpha_new;
    %lambda_old = lambda;
    
    % save current iterate
    iterates(:,:,it) = real( ifdct( wOld, ccindex, N, N ) )';
    
    if( exit_flag )
        return;
    end
end

%B = real( ifdct( wNew, ccindex, N, N ) );

if( printLog )
    %fprintf(1,'it=%5d;\terr=%e;\tfval=%e;\tl1-norm=%e;\tstep=%e;\tdiff=%e;\n',it,err(it),fval,l1norm,alpha,diff);
    fprintf(1,'Eof Main Computation\n\n');
end
%Eof

% =========================================================================
function [Y,angles,truncate,imgSize,maxit,intialGuessStr,constThreshValue] = parseinputs( varargin )

narginchk(2,7);

if( nargin>=1 ) 
    Y = varargin{1};
end

if( nargin>=2 ) 
    angles = varargin{2};
end

truncate = false;
if( nargin>=3 ) 
    truncate = varargin{3};
end

imgSize = [];
if( nargin>=4 ) 
    imgSize = varargin{4};
end

maxit = 50;
if( nargin>=5 ) 
    maxit = varargin{5};
end
    
intialGuessStr = 'zero';
if( nargin>=6 )
    if( ~ischar(varargin{6}) )
        error('Initial guess parameter must be a string!');
    end
    intialGuessStr = varargin{6};
end

constThreshValue = 1;
if( nargin>=7 )
    constThreshValue = varargin{7};
end

% =========================================================================
function keylistener(src,eventdata)

global exit_flag

if strcmp(eventdata.Key,'x')
    exit_flag = true;
    disp('Exit Main routine');
end
